/*  sdftrinf.cpp  by K.Tsuru  */
// function ID 3201 DRADIX
/*************************************************************************
SDouble class
It provides the method to evaluate sin(x), cos(x) or tan(x) functions.
It returns "sign" which must be attached. When sign=0 "y" has an exact
value (-1.0, 0, or 1.0). If you set *func=COS_CALC, SIN_CALC or TAN_CALC,
the function to be used can be gotten by enum number
{COS_CALC, SIN_CALC, TAN_CALC, COT_CALC}.
Refering the value of "*func" the value of sin(y), cos(y), tan(y) or cot(y)
is calculated and if the returned value "sign" is negative you have to change
the sign.
sign == 0 : sin(x) or cos(x) = y ----- y = 1.0 , 0 or -1.0, [+|-]sqrt(0.5)
sign != 0 : sin(x) or cos(x) = sign*(*func)(y)  (|y|<pi/4 )

add "quadrant" of "x" Ver. 2.18 
***************************************************************************/
#ifndef SN_H
#include "sn.h"
#endif

int GetTriCalcMethod(const SDouble& x, SDouble& y, int* func, int* quadrant) {
    const SDouble pi2  = MPi2();        // pi/2
    const SDouble pi4 = MPi4();     // pi/4
    const SDouble m2pi = M2Pi();        // 2/pi
    if(x.Sign(3201) == 0){
        if(quadrant != NULL) *quadrant = 0;
        y = (*func == COS_CALC) ? 1.0 : 0.0;
        return 0;
    }
    double dblX = doubleD(x, 0);
    if(fabs(dblX) < M_PI_4){  // |x| < pi/4
        if(quadrant != NULL) *quadrant = dblX > 0 ? 1 : 4;
        y = Dabs(x);        return (*func == COS_CALC) ? 1.0 : x.Sign(); // cos(x) >0, sin(x) has the same sign as x
    }
/*
It reduces |x| < pi/4.
When |x| has a large value it raises up the precision by the figures of exponent.
It will be rare that |x|>=DRADIX, the recalculation of pi is necessary, then
when PreferSpeed is ON it does not raise the precision.
*/
#ifndef NDEBUG
    assert(x.NetRdxExp() >= 0);
#endif
    int up = 2;
    uint upPrec = (uint)max(up, x.NetRdxExp()); // x.NetRdxExp() ~ 10 , |x| is very large
    if(SNManager::PreferSpeed() == ON) upPrec = 0;

    RealSize C;
    const SDouble HALF(0.5);
    SDouble n, Y;
    if(upPrec) C.SetEffFig(x.EffFig()+upPrec, C.TEMP_EXTEND);   //temporally raise up the precision

    Y = Modf(x*m2pi, n);// |x|/(pi/2) = n+Y (n:integer, 0<=|Y|< 1.0)  // very fast
    if(Y.IsPossibleToRound(2)) Y.Round(); 

    if(quadrant != NULL) {  //Ver. 2.18 (May 9, 2007)
        SLong N;
        N = n; 
        int r4 = N(0) % 4, qd = 0;  // r4 = |n| % 4
        
        if( n.Sign() >= 0 && Y.Sign() > 0) qd = r4 + 1;     // n%4 + 1
        else if( n.Sign() <= 0 && Y.Sign() < 0) qd = 4 - r4;// 4 - |n|%4
        *quadrant = qd; // <= 4
    }
    int c = DDCompare(Y, HALF); //It comapres the absolute values.
    if(c > 0){  // |Y| > 0.5, c> 0
        if(Y.Sign() > 0) n += ONE;  // Y >  0.5, n+Y = (n+1) + (Y-1)
        else n -= ONE;              // Y < -0.5, n+Y = (n-1) + (Y+1)
    } else if(c==0) {// Y=1/2 |x| = pi/4
        y = Sqrt(0.5);
        if((x.Sign() < 0) && (*func == SIN_CALC)) y.ChangeSign();
        return 0;
    }
    y = x - n*pi2;      // y = x - n*pi/2;, |y| <= pi/4 y = pi/4 - pi/2 = -pi/4
    if(upPrec){     C.SetEffFig(0);     y.Reform(3201);     }
#ifndef NDEBUG
 assert(n.RdxExp() >= 0);
#endif
//It determines the sign and function to be used. "n" is an integer.
    int sgn, uf = *func, nsgn = n.Sign(3201);
    uint n1pos = (uint)n.RdxExp();
    fType nL = n(n1pos), k; //the first position of n
                            //The operator() returns zero for the argument without range.
    // x = n*pi/2 + y
    // sin(x) = cos(n*pi/2)*sin(y) + sin(n*pi/2)*cos(y)
    // cos(x) = cos(n*pi/2)*cos(y) - sin(n*pi/2)*sin(y)
    if(nL & 1){ // n = 2k+1, x = (2k+1)*pi/2 + y
        k = (nL -1)/2;
        if(uf == COS_CALC){
            // cos(x) = - sin(n*pi/2)*sin(y)
            *func = SIN_CALC;   // cos(x) = (-1)^(k+1)sin(y)
            sgn = (k & 1) ? nsgn : -nsgn;
        } else if(uf == SIN_CALC){
            // sin(x) = sin(n*pi/2)*cos(y)
            *func = COS_CALC;   // sin(x) = (-1)^k*cos(y)
            sgn = (k & 1) ? -nsgn : nsgn;
        } else {                // tan(x) = -cot(y)
            *func = COT_CALC;   sgn = -1;
        }
    } else {    // n = 2k, x = k*pi + y
        k = nL/2;
        if(uf == COS_CALC){
            // cos(x) = cos(k*pi)*cos(y)
            sgn = (k & 1) ? -1 : 1; // cos(x) = (-1)^k*cos(y)
        } else if(uf == SIN_CALC){
            // sin(x) = cos(k*pi/2)*sin(y)
            sgn = (k & 1) ? -1 : 1; // sin(x) = (-1)^k*sin(y)
        } else sgn = 1;             // tan(x) = tan(y)
    }
/*
When |x| is very large and very close to n*(pi/2), in the evaluation y=x-n*(pi/2)
a cancelation will occure, then "y.IsMLT(one)" must not be used.
*/
    if(y.IsMLT(x)){ // y = 0, x = n*(pi/2), exact value. It does not use "y.IsMLT(one)".
        if(*func == COS_CALC)      y = sgn;
        else if(*func == COT_CALC) y.SetError(y.DOMAIN_ERR,"Tan",3201); // 1/tan(n*pi) = 1/0
        else y.SetZero();   // y = 0 for sin(y)
        sgn = 0;
    }
    return sgn;
}